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ABSTRACT 

Direct collapse of supermassive stars (SMSs) is a possible pathway for generating supermassive black 
holes in the early universe. It is expected that an SMS could form via very rapid mass accretion with 
M* ~ 0.1 — 1 Mq yr~^ during the gravitational collapse of an atomic-cooling primordial gas cloud. 
In this paper we study how stars would evolve under such extreme rapid mass accretion, focusing 
on the early evolution until the stellar mass reaches 10"^ Mq. To this end we numerically calculate 
the detailed interior structure of accreting stars with primordial element abundances. Our results 
show that for accretion rates higher than 10~^ Mq yr~^, stellar evolution is qualitatively different 
from that expected at lower rates. While accreting at these high rates the star always has a radius 
exceeding 100 Rq, which increases monotonically with the stellar mass. The mass-radius relation for 

1/2 

stellar masses exceeding ^100 Mq follows the same track with i?* oc M, in all cases with accretion 
rates > lO'^ Mq yr"!; at a stellar mass of 10^ Mq the radius is ~ 7000 Rq (~ 30 AU). With higher 
accretion rates the onset of hydrogen burning is shifted towards higher stellar masses. In particular, 
for accretion rates exceeding Af* ^0.1 Mq yr"-'^, there is no significant hydrogen burning even after 
10^ Mq have accreted onto the protostar. Such "supergiant" protostars have effective temperatures 
as low as Tefj — 5000 K throughout their evolution and because they hardly emit ionizing photons, 
they do not create an HII region or significantly heat their immediate surroundings. Thus, radiative 
feedback is unable to hinder the growth of rapidly accreting stars to masses in excess of 10'^ M©, as 
long as material is accreted at rates M* > 10~^ Mq yr~^. 

Subject headings: cosmology: theory - early universe - galaxies: formation - stars: formation - 
accretion 
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1. INTRODUCTION 

Recent observations reveal that supermassive black 
holes (SMBHs) exceeding 10^ Mq already existed in the 
universe less than 1 Gyr after the Big Bang (e.g., iFanI 
[2006tlMortlock et al.l[201lLlTreister erall[20ll . The ori- 
gins of such SMBHs must be intimately related to struc- 
ture formation in the early universe. Some scenarios on 
the birth and growth of SMBHs postulate the existence of 
remnant B Hs from Population H I (Pop HI) stars as thei r 
seeds (e.g.. iMadau fc Reesll^OOH: [Schneider et al.ll2002[ ). 
For several decades theoretical studies have predicted 
that the majority of Po p III stars were very n iassive, 
exceeding 100 Mq fe.g.. iBromm fc LarsonI [20041 ). Pop 
HI stars more massive than 300 Mq end their lives by 
direc tly collapsing to form BHs fe.g.. iHeger fc WooslevI 
120021) . If such a ~ 100 M© BH grows via continuous 
mass accretion at the Eddington limited rate, its mass 
barely attains 10^ Mq in 1 Gyr. 

This scenario, however, has recently been challenged. 
First, it is suspected that most Pop HI stars were much 
less massive than previously thought. A circumstellar 
disk forming after the cloud's collapse easily fragments 
due to gravitational i nstability and c ould pr oduce mul- 
tiple protostars (e.g.. iMachida et al. 2008; Stacv et al.l 
120101 Iciark et all 120111 ). The final stellar masses would 
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be reduced as the accreting gas is shared by multiple 
stars. Moreover, strong stellar UV light creates an HII 
region around the protostar when the stellar mass ex- 
ceeds a few xlO Mq. The resulting feedback terminates 
the growth of Pop HI protostars via mass accretion at a 



few xlO MQ_(e 



:S^ 



l2011bt[Stacvet al.l 



McKee fc Tal3 [20081: iHosokawa et al.l 
20121) . A large amount of gas would 



be expelled from the dark halo due to the expansion 
of HII regi ons and the onset o f core-collapse super- 
novae (e.g.. iWhalen et a l. 2004': iKitavama etall 120041 : 
iKitavama fc Yoshidal 1200 5). which quenches the supply 
of gas to any remnant BH. Even if a BH gets some 
gas supply, radiative feedback from the BH accretion 
disk could regulate mass accre t ion onto the BH-disk sys- 
tem fe.g.. [ Alvarez et all 120091 : iMilosavlievic et~all 120091 : 
iJeon et al.ll2'oTlf r 

Another pathway for generating SMBHs is BH binary 
mergers. However, this process could be also limited due 
to the strong recoil resulting from g: r avitational wave 
emission (e.g., iCampanelli et al.l 120071 : iHerrmann et al.l 
|2007). It is not straightforward that a seed BH < 100 Mq 
can grow to a ~ 10^ Mq SMBH within 1 Gyr of its birth. 

An alternative possibility is that massive BHs exceed- 
ing 10^ Mq form d irectly in some rare o ccasions in the 
primordial gas (e.g. . IBromm fc Loebll2003l ). The primary 
cooling process in the primordial gas is line emission of 
molecular hydrogen. However, the thermal evolution of a 
gravitationally collapsing cloud can change significantly, 
if this cooling process is suppressed, for example, due to 
photodi ssociation of n i olecules by stro ng back ground ra- 
diation (|Omukaill200lt lOh fc Haimari .2002: Shane et al.l 
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I2010tllnavoshi fc Oniukai2011f) or coUisiona l dissociation 
in dense shocks (|lnavoshi fc Oinuka]||2012[ ). If the dark 
halo is sufficiently massive (> 10^ Mq), the baryonic gas 
can contract with atomic hydrogen cooling even with- 
out molecular hydrogen. The collapse proceeds nearly 
isothermally at ~ 8000 K. Without efficient molecular 
cooling fragmentation is suppressed and single or binary 
protostars form within dense cloud cores of > 1 0^ Mq 
ijBromm fc Loebl [20031: IRegan fc HaehnelH [2009I ). The 
protostar's mass is initially ^ 10~^ Mq but quickly in- 
creases via mass accretion. The expected accretion rates 
are 0.1 — 1 M© yr~^, more than 100 times higher than 
the standard value ~ 10""^ M© yr~^ expected for Pop III 
star formation. The stellar mass could reach 10^~^ M© 
in ~ 1 Myr with this very rapid mass accretion. General 
relativity predicts that such a supermassive s tar (SMS) 
becomes unstable (e.g. iChandrasekhail 119641 ) and col- 
lapses to form a BH, which subsequen tly swallows most of 
the su rrounding stellar material fe.g.. lShibata fc Shapiro! 
I2002f) . Some authors are exploring a different picture, 
whereby only a central part of the SMS collapses to 
form a ^ 100 Af© BH and heat input from the accreting 
BH inflates the outer e nvelope of the SMS ("quasi-star" 
Begelman eTall [2001 [200l iBegelmanl [2010t iBall et al.l 



20111: iDotan et al.ll201lD ~ 



However, we only have limited knowledge on how 
stars evolve un der such extreme conditions of rapid 
mass accretion. IBegelmanl (j2010( ) predicts that, based 
on simple analytic arguments, such stars have a very 
different structure from their main-sequence counter- 
parts. Stellar evolution at lower accretion rates 
M* ^ 10"^ Mq yr~^ has been studied in de- 
tail by nu merically solv ing the stellar interior struc 



ture (e.g.. IStahler et all Il9 86: Omukai fc Palla 
20031 lOhkubo et al.l pOOOTlHosokawa fc Omukai 
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Hosokawa et al.l I201ol) . lOmukai fc Fallal (|2001L 12003) 
showed that rapid mass accretion with Af, > 4 x 
10~^ Mq yr~^ causes the protostar's abrupt expan- 
sion before its arrival to the Zero-Age Main Sequence 
(ZAMS). Further comprehensive studies on stellar evo- 
lution with rapid mass accretion are indispensable for 
considering their radiative feedback and observational 
signatures (e.g., Johnson et al- 2011). 

We present here our first results of this sort, whereby 
as a first step, we study the early evolution up to a stellar 
mass of 10'^ Af©. Our results show that rapid accretion 
with M^, > 10^^ Af© yr~^ causes the star to bloats up 
like a red giant. The stellar radius increases monotoni- 
cally with stellar mass and reaches ~ 7000 Rq (~ 30 AU) 



tion rates previously studied, the mass-radius relation in 
this phase is almost independent of the assumed accre- 
tion rate. Such massive "super-giant" protostars could 
be the progenitors that eventually evolve to the observed 
SMBHs in the early universe. 

The organization of this paper is as follows. First, 
we briefly review our numerical method and summarize 
the calculated cases in Section [2l We describe our re- 
sults in Section |3l we flrst focus on the fiducial case with 
Af = 0.1 Af© yr~^ and then examine effects of varying 
accretion rates and boundary conditions. Finally, sum- 
mary and discussions are described in Section |4l 



2. NUMERICAL MODELING OF ACCRETING STARS 
2.1. Method 

We calculate stellar evolution with mass accretion us- 
ing the numerical codes developed in our previous work 
(see [OiT Lukai & Palla 2003; Hosokawa fc Omukai 200^ 
iHosokawa et al.ll2010l for details) . The four stellar struc- 
ture equations, i.e., equations of continuity, hydrostatic 
equilibrium, energy conservation, and energy transfer, 
including effects of mass accretion are solved assuming 
spherical symmetry. We focus on the early evolution un- 
til slightly after the ignition of hydrogen fusion in this 
paper. To this end, the appropriate nuclear network for 
the thermo-nuclear burning of deuterium, hydrogen, and 
helium is considered. 

The codes are designed to handle two different outer 
boundary conditions for stellar models: shock and pho- 
tospheric boundaries. The shock boundary condition 
presupposes spherically symmetric accretion onto a pro- 
tostar, whereby the infiow directly hits the stellar sur- 
face and forms a shock fr ont (e.z.. IStahler et all 119801 : 
IHosokawa fc Omukail 120091) . We solve for the structure 
of both the stellar interior and outer accretion flow. 
With this boundary condition the photosphere is lo- 
cated outside the stellar surface where the accretion flow 
is optically thick to the stellar radiation. The photo- 
spheric boundary condition, on the other hand, presup- 
poses a limiting case of mass accretion via a circum- 
stellar disk, whereby accretion columns connecting the 
star and disk are geometrically co mpact and most of the 
stellar surface radiate s freely (e.g.. lPalla fc Stahlerill992l : 
IHosokawa et al.l 120101 ) . In this case we only solve the 
stellar interior structure without considering details of 
the accretion flow; the location of photosphere always 
coincides with the stellar surface. 

The different outer boundary conditions correspond to 
the two extremes of accretion flow geometries, or more 
specifically to different thermal efficiencies of mass accre- 
tion, which det ermine the specific ent ropy of accreting 
materials (e.s... [Hosokawa et al.|[2QllaD . The accretion 
thermal efficiency controls the entropy content of the 
star, which determines the stellar structure. With the 
shock boundary condition, the accreting gas obtains a 
fraction of the entropy generated behind the shock front 
at the stellar surface. The resulting thermal efficiency is 
relatively high ( "hot" accretion) . 

For the photospheric boundary condition, on the other 
hand, the accreting gas is assumed to have the same en- 
tropy as in the stellar atmosphere. The underlying idea 
is that, when the accreting gas slowly approaches the star 
via angular momentum transport in the disk, its entropy 
should be regulated to the atmospheric value. This is 
a limiting case of thermally inefficient accretion ( "cold" 
accretion). In general, with even a small amount of angu- 
lar momentum, mass accretion onto the protostar would 
be via a circumstellar disk, perhaps with geometrically 
narrow accretion columns connecting the disk with the 
star. For extremely rapid mass accretion, however, the 
innermost part of the disk becomes hot and entropy gen- 
erate d within the dis k is advected into the stellar interior 
(e.g.. [Popham et al.i 1993) . Thus, cold accretion as en- 
visioned for the photospheric boundary condition is not 
appropriate fo r the case of rapid mass accretion (also se e 
discussions in lHartmann et aL[|1997l: [Smith et al.l 120 111) . 
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TABLE 1 

Cases Considered 



Case 



M. {Mq yr-l) M.,0 (M©) i?,,o {Rq) 



Notes 



References 



MDleO 


1.0 


MD3el 


0.3 


MD3el-HC10 


0.3 


MD3el-HC50 


0.3 


MDlel 


0.1 


MD6e2 


0.06 


MD3e2 


0.03 


MD6e3 


0.006 


MDle3 


0.001 



2.5 

2.0 
2.0 (10) 
2.0 (50) 

2.0 

1.0 

1.0 

1.0 

0.05 



298.4 

238.3 
238.3 (437.1) 
238.3 (826.8) 

177.8 

118.6 

95.2 

52.3 

12.5 



shock 
shock 



photo. BC at M« 

photo. BC at M* 

fiducial case 



10 Mq 

50 Mq 



also see Omuka^&^alla (2003) 
also see ^muka^_&_^alla (2003) 



Sec.O 

Sec. O EH 

Sec. [331 

Sec. [331 

Sec. [3ll] 1321 

Sec. (321 

Sec. [3:21 

Sec. (321 

Sec. (32] 13:21 



Note. — Col. 2: mass accretion rate, Col. 3: initial stellar mass, Col. 4: initial radius. For cases MD3el-HC10 and 
MD3el-HD50 the values when the boundary condition is switched is given in parentheses. 



We therefore expect that the shock boundary condition 
is a good approximation for the extremely high accre- 
tion rates considered here and mostly focus on stellar 
evolution with the shock boundary condition. We also 
consider a few cases with the photospheric boundary con- 
dition for comparison to test potential effects of reducing 
the accretion thermal efficiency (also see Sec. l2.2l below). 

2.2. Cases Considered 

The cases considered are summarized in Table 1. In 
this paper, we only consider the evolution with constant 
accretion rates for simplicity. The adopted accretion 
rates range from 10~^ M© yr~^ to 1 Mq yr~^. Stellar 
evolution with accretion rates less than 10~^ Mq yr~^ 
(cases MDleS and MD6e3) ha s been studied in de- 
tail in our pre vious work (e.g., lOmukai fc Pallal 120031 : 
iHosokawa fc O mukai 2000). As described in Section [23] 
above, we calculate the protostcllar evolution assum- 
ing shock outer boundary conditions for most of the 
cases. Cases MDBel-HClO and MD3el-HC50 are the 
only exceptions, whereby we switch to the photospheric 
boundary condition after the stellar mass reaches 10 Mq 
and 50 Mq, respectively, at a constant accretion rate 
0.3 Mq yr~^. The underlying idea for switching the 
boundary at some moment is that the specific angu- 
lar momentum of the inflow and thus the circumstellar 
disk grows with time, which reduces entropy of the ac- 
creted matter. The higher mass at switching point corre- 
sponds to the higher angular momentum of the parental 
core. As in our previous work, we start the calcula- 
tions with initial stellar models constructed assuming 
that the stellar interior is in radiative equilibrium (e.g., 
IHosokawa fc Omukaill2009l ) . We adopt a slightly higher 
initial stellar mass of ~ IM© for stability reasons. The 
calculated initial stellar radii are also summarized in Ta- 
ble [H 



3. RESULTS 



3.1. 



Evolution in the Fiducial Case 

(Ah = 0.1 Mq yr^i; 

We first consider the fiducial case (MDlel), whereby 
the stellar mass increases with the constant accretion rate 
of M* = 0.1 Mq yr~^. The calculated evolution of the 
stellar interior structure is presented in Figure [TJ We see 
that the stellar radius is very large and increases mono- 
tonically with the stellar mass. The stellar radius exceeds 
10^ Rq when the stellar mass is M„ ~ 45 Mq, and reaches 
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Fig. 1. — Evolution of the stellar interior structure for the fiducial 
case, whereby the stellar mass increases at a rate of 0.1 Mq yr~^ 
(case MDlel). The thick solid line depicts the stellar surface, which 
is the position of the accretion shock front. The dotted lines show 
the radial positions of the mass coordinates of M = 3, 10, 30, 
100, and 300 Mq. The dot-solid line indicates the radial position 
within which 70 % of the stellar mass is enclosed. The gray-shaded 
areas represent the convective layers. The hatched areas indicate 
layers of active nuclear burning, where the energy production rate 
exceeds 10% of the steady rate 0.1LD,at/A/* for deuterium burning 
(see eq. IIOII . and 0.1L*/M» for hydrogen burning. The blue solid 
line shows the evolution of the radius of a star accreting material 
at 10"'^ Mq yr~^ for comparison. 

~ 6500 Rq at M* ~ 10^ Mq. This evolution differs qual- 
itatively from that calculated assuming a lower accretion 
rate 10""^ Mq yr~ ^ as depicted in Figur e [1] by the blue 
fine (see also e.g.. lOmukai &: Pallal 120031 ). At the lower 
accretion rate the stellar radius initially increases with 
stellar mass but begins to decrease at M* > 6 Mq. 

A key quantity for understanding the contrast between 
the two cases is the balance between the two character- 
istic timescales: the accretion timescale 



t. 



M, 
AL 



and Kelvin-Helmholtz (KH) timescale 

GM^ 



R.tL„ 



(1) 



(2) 



where _R* and i* are the stellar radius and luminosity, 
and G is the gravitational constant (e.g., .Stabler et ahl 
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Fig. 2. — Evolution of the stellar surface luminosity L» (dashed 
line) and maximum luminosity within the star Lmax (solid line) 
for the cases with 10~^ Mq yr~^ (MDleS; upper panel) and 
0.1 Mq yr~^ (fiducial case MDlel; lower panel). The mass- 
luminosity relations given by equations ((Sj, lO and JSj are shown 
with the thin green, blue, and red lines, respectively. In each panel 
the vertical dot-dashed line (magenta) indicates the epoch when 
the accretion time is equal to the KH time. 



\T9M\(h^^^i fc Pallal[2nn3l : IHosokawa fc OmukaillMifll ). 
In the early stage, during which the steUar radius in- 
creases with mass, the timescale balance is iacc ^ ^kh 
and radiative loss of the stellar energy is negligible ( adi- 
abatic accretion stage). However, the radiative energy 
loss becomes more efficient as the stellar mass increases. 
This is because opacity in the stellar interior, which is 
due to the free-free absorption following Kramers' law 
K ex pT~^-^, decreases and the stellar luminosity L* in- 
creases with the stellar interior temperature (and thus 
with its mass M^). The upper panel of Figure [5] indeed 
shows that the maximum luminosity within the star L^ax 
increases as a power-law function of M^ . 

This increase of L, is consistent with the analytic scal- 
ing relation for radiative stars with Kramers' opacity, 

L ex Af"/^i?::^/^ (e.g- IHavashiet al.l[r962) . Our nu- 
merical results are well fitted by the analytic relation 



Ln 



0.6 L, 







Ma- 



11/2 



Rq 



-1/2 



(3) 



The increase of L^ causes an inversion of the timescale 
balance to iacc > ^kh at low accretion rates. The star 
contracts by losing its energy (KH contraction stage), 
which is seen for Af* > 6 Mq. The opacity in the stellar 
interior has fallen down to the constant value of electron 
scattering. Figure [2] shows that, in this stage, luminos- 
ity takes its maximum value at the stellar surface and 
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Fig. 3. — Evolution of several physical quantities for the fiducial 
case with 0.1 Mq yr~^ (MDlel). Top panel: Comparison be- 
tween the accretion timescale tacc (dashed line) and KH timescale 
txH (dotted line). The vertical magenta dot-dashed line indicates 
the epoch when iKH is equal to tacc The thin solid line repre- 
sents 100 times of the stellar free-fall timescale tff = ^37r/32Gp, 
where p is the average mass density of the star. The fact that tff is 
much shorter than (kh and tacc verifies the hydrostatic balance as- 
sumption implicit in the stellar structure equations. Middle panel: 
Evolution of the accretion luminosity Lace (dashed line), stellar 
luminosity L* (dotted line), and total luminosity Ltot = iacc -f L« 
(solid line). The red line indicates the Eddington luminosity at 
each stellar mass. Bottom panel: Evolution of the radial posi- 
tions of the photosphere ijpij(ij0) (dashed line) and stellar surface 
R*{Rq) (black solid line). The red solid and dashed lines denote 
the analytic formulae for these radii by Stahlcr ct al. (1986) feas[6l 
andO. The evolution of the stellar effective temperature T^ff (X) 
is also shown with the blue line. 



increases as 



oc 



ML 



which is valid for the 



constant opacity cases (e.g.. iHavashi et aLlll962[ ). The 
relation 

3 



Lr, 



10 L 







M^ 



(4) 



roughly agrees with our results. Temperature at the stel- 
lar center increases during the KH contraction stage. 
Hydrogen burning finally begins and the stellar radius 
begins to increase following the mass-radius relation of 
ZAMS stars for M, > 50 Mq. Figure [2] shows that the 
stellar luminosity gradually approaches to the Eddington 
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limit 



iEdd(M*) ~ 3.8 X 10^ ie 



Ah 



100 M, 



(5) 



By contrast, there is no contraction stage for the case 
with a much higher accretion rate M* =0.1 Mq yr~^ 
(MDlel). Nevertheless, the evolution of the timescales 
still follows the picture described above (Fig. [3] a). We 
see that the timescale balance changes from tacc < ^kh 
to iacc > ^KH at M* ~ 40 Mq. The protostar is in the 
adiabatic accretion stage for M, < 40 Mq. The accre- 
tion luminosity Lace = GM^M^/R^ at this stage is much 
higher than the stellar luminosity L*, since the luminos- 
ity ratio La^cc/ L* is equal to the tim escale ratio tKn/iacc 
by definition. iStahler et al.l (|1986[ ) derived the approxi- 
mate analytic formulae describing radial positions of the 
stellar surface i?, and photosphere i?ph (located within 
the accretion flow) during the adiabatic stage: 



i?, ~ 26 R 







Ma- 



0.27 



M, 



0.41 



10-3 Mq yr-1 



(6) 



i?ph ^ 1.4i?,. (7) 

The bottom panel of Figure |3] shows that these 
formulae still agree with our numerical results with 
0.1 Mq yr^^. Equation © shows that the stellar ra- 
dius is larger for the higher accretion rate at the same 
mass. The larger radius is due to the higher specific 
entropy of accreting material and the resulting higher 
entro py content of the star (e.g., iHosokawa fc Omukail 
|2009( ). Comparing two stars of the same mass, the one 
with a larger radius has a lower interior temperature, 
which then implies a higher opacity due to the strong T- 
dependence of Kramers' law k cx pT~^'^ . For this reason 
the adiabatic accretion stage is prolonged up to higher 
stellar masses for higher accretion rates. The stellar lu- 
minosity L* thus increases with stellar mass. Figure [5] 
(b) shows that the evolution of the stellar maximum lu- 
minosity Lmax still obeys equations ([3]) - ([5]). Unlike in 
the case for 10"^ Mq yr~^, however, it is only after imax 
approaches the relation of Lmax k M^ (eq. H]) that ^kh 
becomes equal to iacc- The rapid heat input via mass 
accretion prevents the star from losing internal energy 
until the stellar luminosity becomes sufficiently high. 

The fact that ^kh is shorter than iacc for M, > 40 Mq 
indicates that most of the stellar interior is contract- 
ing, as shown by the trajectories of the mass coordinates 
(dashed lines in Figure [T|). The figure also shows that the 
bloated surface layer occupies only a small fraction of the 
total stellar mass. When the stellar mass is ~ 300 Mq, 
for example, the layer which has 30 % of the total mass 
measured from the surface covers more than 98 % of the 
radial extent. The star has a radiative core surrounded 
by an outer convective layer. Although the convective 
layer covers a large fraction of the stellar radius, even 
the central radiative core is much larger than a ZAMS 
star with the same mass (compare with the blue curve 
for Af, > 40 Mq in Fig. [T]). Figure H shows the radial 
distributions of physical quantities, i.e., specific entropy, 
luminosity, temperature, and density, in the stellar inte- 
rior. We see that the specific entropy is at its maximum 
value near the boundary between the radiative core and 



convective layer. The stellar entropy distribution is con- 
trolled by the energy equation. 



T 



M 



f dL 
\dM,, 



(8) 



where s is the specific entropy and e is the energy produc- 
tion rate by nuclear fusion. For Af* > 40 Mq most of the 
stellar luminosity comes from the release of gravitational 
energy. In fact, as seen in Figure U (b), {dL/dM)t > 
in the radiative core, which means that the internal en- 
ergy of the gas is decreasing. The local luminosity in the 
radiative core is close to the Eddington value given by 
equation ([5]), using the mass coordinate M rather than 
stellar mass M* (Fig. |3]b). Note that the opacity in the 
radiative core is only slightly higher than that expected 
from electron-scattering alone. 

In the outer parts of the star, where temperature and 
density are lower, however, opacity is higher than in the 
core because of bound- free absorption of H, He atoms and 
the H^ ion. Energy transport via radiation is inefficient 
there, and a part of the energy coming from the core is 
carried outward via convection. Figure U shows that the 
surface convective layer lies in the temperature and den- 
sity range where T < 10^ K and p < 10"^ cm^^^ which is 
almost independent of the stellar mass. We also see that 
the specific entropy is not constant over the convective 
layer, decreasing toward the stellar surface (a so-called 
super- adiabatic layer). This is because convective heat 
transport is inefficient near the stellar surface. A part of 
the outflowing energy is absorbed there, as indicated by 
the fact that the surface layer has a negative luminosity 
gradient {dL/dM)t < 0. This explains the high specific 
entropy in the outer convective layer. 

The outermost part of the star has a density inversion, 
i.e. the density increases toward the stellar surface. Here, 
opacity assumes very high values because of H~ absorp- 
tion. Radiation pressure is so strong that the hydrostatic 
balance is not achieved only with gravity; the additional 
inward force by the negative gas pressure gradient helps 
maintain the hydrostatic structure. Note, however, that 
this density inversio n could be unstable in r ealistic multi- 
dimensions (see e.g.. iBegelman et "aLll2008[ ). 

Although deuterium burning is ignited when the stellar 
mass is ~ 50 Mq, its influence on the subsequent evolu- 
tion is negligible (Fig. [T]). The total energy production 
rate by deuterium burning is approximately 



= 1.5 X 10^ Lq 



(9) 



M, 



[D/H] 
0.1 Mq yr-1 / V2-5 X 10-5 



where djj is the energy available from deuterium burn- 
ing per unit gas mass. Since this is much lower than 
the Eddington luminosity iEdd in the mass range con- 
sidered, energy production by deuterium burning con- 
tributes only slightly to the luminosity in the stellar inte- 
rior. Hosokawa & Omukai (2009) showed that deuterium 
burning influences the stellar evolution only when the ac- 
cretion rate is low M* < 10"^ Mq yr"!. 

Figure 131(c) shows that temperature in the stellar inte- 
rior increases with total mass. The central temperature 
reaches 10* K when the stellar mass is ~ 600 Mq. Soon 
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Fig. 4. — Radial profiles of the specific entropy (panel a), lumi- 
nosity (6), temperature (c) and gas mass density (d) in the stellar 
interior for the fiducial case with 0.1 Mq yr~^ (MDlel). The 
profiles when the stellar mass is 200 Mq, 600 Mq, and 10^ Mq 
are shown in each panel. The magenta parts indicate the lay- 
ers are convective. The thin red line in panel (b) represents the 
Eddington-limit luminosity as a function of the mass coordinate 
M. 



after that, hydrogen burning begins and a central con- 
vective core develops (Fig. [T]). This convective core can 
also be seen in the radial profiles for the 10"^ Mq model 
in Figure 21 (indicated by magenta). The luminosity pro- 
file tells that most of the energy produced by hydrogen 
burning is absorbed within the convective core. The star 
still shines largely by releasing its gravitational energy 
even after hydrogen ignition. 

When the stellar radius is sufficiently large, the ac- 
cretion flow reaches the stellar surface before becoming 
opaque to the outgoing stellar light. In fact, soon after 
the end of the adiabatic accretion stage, the accreting en- 
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Fig. 5. — Evolution of the protostellar radius for various ac- 
cretion rates. Upper panel: the different curves represent the 
cases with M, = 10~^ Mq yr'^ (case MDleS, black), 6 X 
10-3 Mq yr-l (MD6e3, blue), 3 X lO'^ Mq yr"! (MD3e2, red), 
and 6 X 10~^ Mq yr'^ (MD6e2, magenta). The open and filled 
circles on each curve denote the epoch when (kh = iacc and 
when the central hydrogen burning begins, respectively. Lower 
panel: same as the upper panel but for higher accretion rates of 
6 X 10-2 Mq yr-i (MD6e2, magenta), 0.1 Mq yr-^ (MDlel, red), 
0.3 Mq yr-l (MD3el, blue), and 1 Mq yr"! (MDleO, black). The 
case MDe2 is illustrated in both panels as a reference. For the cases 
with 0.3 Mq yr-i and 1 Mq yv'^ (MD3el and MDleO) hydrogen 
fusion has not ignited by the time the stellar mass reaches 10'^ Mq . 
In the both panels the thi n gr een line represents the mass-radius 
relation given by equation l|12p . 



velope remains optically thin throughout (Fig. [3]c). We 
also see that the stellar effective temperature is almost 
constant at Toff ~ 5000 K during this period. In general, 
the stellar effective temperature never assumes a lower 
value due to the strong temperature-d ependence of H^ 
absorption opacity (e.g.. lHavashilll96ll) . Stars that have 
a compact core and bloated envelope (e.g., red giants) 
commonly have an almost constant effective tempera- 
ture, regardless of their stellar masses. 

3.2. Cases with Different Accretion Rates 

We now investigate how stellar evolution changes with 
the accretion rate. Figure \5\ shows the evolution of 
the stellar radius for several cases, including M* = 
10-3 Mq yr-i (case MDleS) and 0.1 Mq yr'^ (case 
MDlel) explained in Section [01 We see that, with ac- 
cretion rates higher than 6 x 10"^ Mq yr"^, the evo- 
lution becomes similar to that of the fiducial case for 
0.1 Mq yr"^ (MDlel); the stellar radius monotonically 
increases with mass. The protostars undergo adiabatic 
accretion in the early stage. As equation ^ shows, the 
stellar radius is larger for higher accretion rates at a given 
stellar mass, say, at M* = 10 Mq. 
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Fig. 6. — Upper panel: Evolution of radial positions of the stellar 
surface (solid) and photosphere (dotted). The black, blue, and 
magenta curves denote cases with A/* = lO""^ Mq yr"-*^ (case 
MDleS), 6 X 10-3 Mq yr"! (MD6e3), and 1 Mq yr"! (MDleO), 
respectively. The mass-radius relation given by equation II12I I is 
plotted with the thin green line. Lower panel: The stellar effective 
temperature for the same cases as in the upper panel. 
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Fig. 7. — Evolutionary tracks in the HR diagram. The dif- 
ferent colors denote the evolution for different accretion rates, 
10-3 Mq yr-l (case MDleS, black), 6 X IQ-^ Mq yr"! (MD6e3, 
blue), 0.1 Mq yr-i (MDlel, red), 1 Mq yr"! (MDleO, magenta). 
For each case the values at the stellar photosphere (Ltot and T^ff ) 
and at the stellar surface (L« and Toff <, = Lt/A-KcrRl) are plot- 
ted with solid and dashed lines, respectively. The thick dashed 
line represent s the loci of non-accreting ZAM S stars taken from 
IMarieo et al] 1I2OOII ') (M, < 100 Mq) and Hr omm et al] 1I2OOII) 
(M« > 100 Mq). The filled circles and open squares on the lines 
mark the positions for M* = 30 Mq, 100 Mq, 300 Mq, and 
10^ Mq in ascending order. 
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Fig. 8. — Comparisons of the interior structure of 10^ Mq stars 
produced by different accretion rates. Radial profiles of mass {panel 
a), specific entropy (6), and luminosity (c) are shown. In each panel 
the solid and dashed curves represent the cases for 0.3 Mq yr"^ 
(MD3ol), 0.1 Mq yr-i (MDlel), 3x10"^ Mq yr"! (MD3e2), and 
10~3 Mq yr~^ (MDle3). The magenta lines mark convective lay- 
ers in the stellar interior. The thin red line in panel (c) represents 
the Eddington luminosity as a function of mass M. 



Adiabatic accretion occurs up to a higher stehar mass 
when the accretion rate is higher. We derive an analytic 
expression describing this dependence in the following. 
As long as the KH timescale ixH is longer than the ac- 
cretion timescale ^acc i we h ave adiabatic accretion, (e.g., 
iHosokawa fc Omuka i 20091. Figure [2] shows that imax, 
whose evolution is well described by simple analytic ex- 
pressions (eqs. [3] and E]), converges to L, at the end 
of the adiabatic accretion stage. Thus, the stellar mass 
when the adiabatic accretion terminates M*_teq can be 
estimated with the relation tacc — ^kh = G'M^/i?*Linax- 
Eliminating imax and i?* with equations ^ and (|6]), we 
obtain 



M: 



+ ,tcq 



14.9 M, 



AL 



0.26 







10-2 Mq yr 



.-1 



(10) 



We have confirmed that the epoch of the timescale equal- 
ity in our numerical calculation is well described by this 
equation. 

Even after ^kh becomes shorter than iacc, the stel- 
lar radius continues to increase for accretion rates > 
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10^^ Mq yr^^. The variations of radii among cases with 
different accretion rates graduahy disappear. The stel- 
lar radii finally converge to a unique mass-radius relation 
with i?, oc Af, ' in all the cases. We can derive the ap- 
proximate mass-radius relation from the following simple 
argument. First, the stellar luminosity is generally writ- 
ten as 

L, = AuRlaT^s. (11) 

where a is the Stefan-Boltzman constant. As Figure |3] 
indicates, the stellar luminosity approaches the Edding- 
ton value iEdd(-M*) for M, > 100 Mq. As explained in 
Section [01 the stellar effective temperature stays at the 
constant value Tcff ~ 5000 K after the inversion of the 
timescales (also see Figs. |6] and [T]). Substituting these 
relations into equation (|lip. we obtain 



i?* 



2.6 X 10-^ R, 



■o 



M, 



100 M, 



1/2 



(12) 



Figure [5] shows that our numerical results approximately 
follow this relation. 

Begelman (2010) also considered stellar evolution with 
very rapid mass accretion using simple analytic argu- 
ments. His model predicts that the stellar radius is pro- 
portional to the mass accretion rate (his eq. 24), which 
does n ot agree with our numerical results. iBegelmanl 
(j2010t ). however, does not take into account the de- 
tailed structure of the outermost layer of the star, where 
H^ opacity is important. The fact that the strong T- 
dependence of H~ opacity keeps the stellar effective tem- 
perature almost constant is essential for our results. 

Cases MD3e2 and MD6e3 model stellar evolution 
at intermediate accretion rates (3 x 10~^ and 6 x 
10^"^ Mq yr^^, respectively) and exhibit a different be- 
havior than the higher accretion-rate cases described 
above. For an accretion rate 3 x 10~^ Mq yr~^ (MD3e2) 
the protostar initially contracts after the adiabatic ac- 
cretion stage. At M* ~ 70 Mq, however, the stellar 
radius sharply increases and ultimately converges to the 
mass-radius relation given by equation ([12]). The case 
with 6 x 10~^ Mq yr"-'^ exhibits an oscillatory behav- 
ior of the stellar radius for M* > 70 Mq. In this 
case the accreting envelope remains optically thick af- 
ter the onset of KH contraction (Figs. |6] and [7|). Its 
photospheric radius still follows the mass-radius relation 

1/2 

i?ph OC M, ' . The effective temperature also assumes 
the constant value Toff ~ 6000 K . These evolutionary 
features for M* < 1Q~^ Mq yr~^ have also been fou nd 
in previous studies fe.g.. lOmukai fc Pallal[20Qlll2003D . 

We have seen that the stellar radius at M* ~ 10'^ Mq 
is almost independent of accretion rate as long as Af» > 
3 X 10~^ Mq yr~^. However, the stellar interior struc- 
ture at this moment is not identical among these cases 
(Fig. [5]) . Although each of these stars has a radiative core 
and a convective envelope, the mass is more strongly cen- 
trally concentrated for the lower accretion rates; the less 
massive envelopes have a higher entropy and inflate even 
more to achieve the same stellar radius. 

The evolution of the stellar maximum temperature is 
helpful for understanding the variation of the stellar inte- 
rior structure (Fig. [9] a) with accretion rate. As equation 
p^ shows, the central part of the star begins to contract 
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Fig. 9. — Evolution of the maximum temperature in the stel- 
lar interior {upper panel) and Eddington ratio itot/iEdd {lower 
panel) with increasing stellar mass. The solid and dashed curves 
alternately represent the cases with different accretion rates, 
10-3 Mq yr-i (case MDleS), 6 X IQ-^ Mq yr"! (MD6e3), 3 X 
10-2 Mq yr-i {MD3e2), 0.1 Mq yr"! (MDlel), and 0.3 Mq yr"! 
(MD3el). 



and release gravitational energy at a lower stellar mass 
for the lower accretion rates. The central temperature 
quickly increases with stellar mass once the KH time be- 
comes shorter than the accretion time. The maximum 
temperature Tmax reaches 10® K at a lower stellar mass 
for the lower accretion rate. After that Tmax assumes an 
almost constant value due to the strong T-dependence of 
the energy production rate of hydrogen burning. For the 
case with 3 x 10~^ M© yr"^ (MD3e2) hydrogen burning 
begins at Af* ~ 200 Mq ; the resulting central convective 
core is seen in the profiles in Figure [S] This feature is 
not seen for the case with 0.3 Mq yr~^, because hydro- 
gen has not yet ignited by the time M* — 10"^ Mq for 
M* > 0.3 Mq yr-i. 

Figure El shows that for M* > 6 x 10"^ Mq yr^^ the 
protos tar cannot reach t he ZA MS stage by KH contrac- 
tion. lOmukai fc Fallal (|2003D pointed out that this is 
because the total luminosity Ltot = i* + iacc becomes 
close to the Eddington limit during the contraction to the 
ZAMS. For the cases with 6 x 10"^ and 3 x 10"^ ^/^ yj.-^ 
(cases MD6e3 and MD3e2), for example, the abrupt ex- 
pansion terminates the KH contraction when the total 
luminosity is ne a rly at the Eddington limit (Fig. IHlb). 
lOmukai fc Fallal (J2003I ) analytically derived the maxi- 
mum accretion rate ~ 4 x 10~^ Mq yr"^ with which 
the protostar can reach the ZAMS following KH con- 
traction. Figure [5] indicates that there is another critical 
accretion rate ~ 6 x 10~^ Mq yr"^, above which the stel- 
lar evolution changes qualitatively; the KH contraction 
stage disappears entirely at higher rates. This critical 
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rate can also be derived froni a similar argument as the 
one above. 

Note that the increase of stellar mass during the 
KH contraction stage is smaller for the case with 3 x 
10~2 Mq yr-i than with 6 x 10"^ Mq yr'^ Extend- 
ing this fact for our critical case, the total luminosity 
would nearly reach the Eddington limit just at the end 
of the adiabatic accretion stage, i.e., when Ikh — iacc- 
Since the opacity in the surface layer is higher than from 
Thomson scattering during this epoch, having the total 
luminosity only slightly lower than the Eddington value 
causes the star to expand. Thus, the condition for the 
critical case is 



10000 



2L„ 



CEdd^Edd, 



(13) 



where Csdd is a factor less than the unity and we have 
used the fact that the total luminosity is written as 
Ltot ^ 2imax when ixH '^ iacc- Using CEdd = 0.25 
as a fiducial value (Fig. [HI b) and equation ^ for L,„ax, 
the stellar mass which satisfies the condition (TTSl) is 



M,,Edd,tcq ^ 21.7 MfT 



CEdd 

0.25 



0.5 



(14) 



On the other hand, equation ((TU)) also gives the stellar 
mass when ^kh — iacc for a given accretion rate. Equat- 
ing M^^tcq and M^.Edd.toq with equations (fTOl) and ([T4|) . 
we obtain the critical mass accretion rate 



M„c.4.7xlO-2Moyr-i(^^) 
which agrees with our numerical results. 



1.9 



(15) 



3.3. Effects of Lower- Entropy Accretion 

We have used the shock outer boundary condition for 
the stellar models presented and discussed above. As 
discussed in Section 12.11 the shock boundary condition, 
which implies that the accreting gas joins the star with 
relatively high entropy, would be valid for cases with the 
very rapid mass accretion considered in this paper. If the 
accred gas had lower entropy, however, the stellar radius 
would be reduced because of the resulting lower entropy 
thoughout the stellar interior. Here, we examine poten- 
tial effects of the colder mass accretion by adopti ng the 
photospheric boundary conditions (e.g., Hosokaw a et al.l 
120101 I2011al) . Figure [10] (a) shows the evolution of the 
stellar radius for three cases with 0.3 Mq yr~^, whereby 
the shock boundary condition is used throughout in one 
case (MD3el), whereas the outer boundary condition is 
changed to the photospheric one for > 10 Mq (MD3el- 
HCmlO) and > 50 Mq (MD3el-HCm50), respectively. 
The stars are still in the adiabatic accretion stage when 
the boundary condition is changed in both cases. The 
different outer boundary conditions do affect the stellar 
evolution. For the case where the photospheric boundary 
condition is adopted at M* = 10 Mq (MD3el-HCnilO), 
for example, the star initially contracts after the bound- 
ary condition is switched at M^, = 10 Mq, and then 
abruptly inflates at M* ~ 45 Mq. The stellar radius 
exceeds 10^ Rq and gradually increases with the stellar 
mass thereafter. In spite of the different behaviors in the 
early stages, the subsequent evolution for M* > 100 Mq 
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Fig. 10. — Effect of reducing the thermal efficiency of mass accre- 
tion {upper panel: stellar radius, lower panel: maximum tempera- 
ture within the star). The same accretion rate of 0.3 Mq yr~ is 
adopted for all three cases presented. The solid line represents the 
evolution with the shock boundary condition, i.e, thermally effi- 
cient or "hot" accretion, throughout (MD3el). The dashed and 
dot-dashed lines show the evolution in cases MD3el-IfC10 and 
MD3el-HC50, where the photospheric boundary condition (i.e., 
thermally inefficient or "cold" accretion) is adopted after the stel- 
lar mass exceeds 10 Mq and 50 Mq, respectively. In the lower 
panel the dot-dashed line is indistinguishable from the solid line. 

is quite similar to that in the case with the shock bound- 
ary condition throughout (MD3el). The evolution when 
the boundary condition switching occurs at M^, = 50 Mq 
(MD3el-IICm50) is much closer to that in the shock- 
boundary case. The uniqueness of the mass-radius re- 
lation for M* > 100 Mq can be explained by the fact 
that the argument leading to the analytic expression, 
equation (I12p . does not assume a specific boundary con- 
dition. When the boundary condition is switched at 
M* = 10 Mq (MD3el-IICml0), the stellar interior tem- 
perature is higher and thus the opacity in the stellar 
interior (oc T^3-5 according to Kramars's law) is lower 
than for the shock-boundary case (MD3el) at the same 
stellar mass. As a result the star begins to release its 
internal energy earlier than for the shock-boundary case. 
Indeed, the timescale equality between ^kh and iacc oc- 
curs at M* ~ 40 Mq , earlier than for the shock-boundary 
case, which occurs at the time of abrupt expansion of the 
stellar radius. 

4. SUMMARY AND DISCUSSIONS 
We have studied the evolution of stars growing via 
very rapid mass accretion with 10~^ Mq yi~^ < M* < 
1 Mq yr~^, which potentially leads to formation of 
SMBHs in the early universe. In contrast to previous at- 
tempts to address this problem, we study the stars' evo- 
lution by numerically solving the stellar structure equa- 
tions including mass accretion. Our calculations show 
that stellar evolution in such cases is qualitatively dif- 
ferent from that expected for normal Pop III star for- 
mation, which proceeds at much lower accretion rates 
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^ 10^'^ — 10^^ Mq yr^^ . Rapid mass accretion causes the 
star to inflate; the stellar radius further increases mono- 
tonically with stellar mass at least up to M* ~ 10^ M©. 
For masses exceeding ~ 100 Mq, the star consists of 
a contracting radiative core and a bloated surface con- 
vective layer. The surface layer, which contains only a 
small fraction of the total stellar mass, fills out most 
of the stellar radius. The evolution of the stellar ra- 
dius in this stage follows a unique mass-radius relation 

i?, ex Afy^ which reaches ~ 7000 i?0(~ 30 AU) at 
M, = 10^ Mq, in ah the cases with > 10"^ Mq yr~^ 
Hydrogen burning begins only after the star becomes 
very massive (M, > 100 Mq)\ its onset is shifted toward 
higher masses for higher accretion rates. With very high 
accretion rates M* > 0.1 Mq yr~^, hydrogen is ignited 
after the stellar mass exceeds 10^ Mq. The stellar radius 

continues to grow as i?* oc M* even after hydrogen ig- 
nition. 

In this paper we have focused on the early evolution 
until the stellar mass reaches 10^ Mq. The subsequent 
evolution remains unexplored because of convergence dif- 
ficulties with the current numerical codes. If the star 
continues to expand following the same mass-radius re- 
lation (|12|) also for A/* > 10^ M©, the stellar radius 
at 10^ Mq would be ~ 400 AU. Since the stellar ef- 
fective temperature remains ~ 5000 K, the star hardly 
emits ionizing photons during accretion. Therefore, it 
is unlikely that stellar growth is limited by the radia- 
tiv e feedback via fo r mation of an HII region a s disc ussed 
by iHosokawa et all ()2011b[ ). I Johnson et "all ()201lD also 
reached an analogous conclusion that UV feedback does 
not hinder SMS formation. In their argument, however, 
the star is assumed to reach the ZAMS and to emit a 
copious amount of ionizing photons, but the expansion 
of the HII region is squelched by rapid spherical infiow. 
They also expected that, as a result of confinement of 
the HII region, strong emission lines reprocessed from 
the ionizing photons (e.g., Ly a and He II) would escape 
from the accretion envelope to be an observational signa- 
ture of these objects. By contrast, our calculations show 
that the stellar UV luminosity and thus the luminosi- 
ties in those lines should b e much weaker than s upposed. 
Note that the argument bv lJohnson et al.l ()2011[ ) assumes 
perfect spherical symmetry, which allows the HII region 
to be confined within the accretion envelope. Given that 
mass accretion will likely occur through a circumstellar 
disk, the HII region should grow toward the polar region 



where the gas dens ity is much lower than the spherical 
accretion flow fe.g. . IHosokawa et al.ll2011b[ ). This should 
be the ca se with t he hig h stellar UV luminosity assumed 
in Johns on et al.l ()2011[ ) . 

Even without stellar radiative feedback, stellar growth 
via mass accretion might be hindered by some other pro- 
cess, e.g., rapid mass loss. Indeed, evolved massive stars 
in the Galaxy (M* '^ 10 — 100 Mq), which have large 
radii (i?* > 100 Rq) and high luminosities close to the 
Eddington limit (L, ~ 10^ Lq), generally have strong 
stellar winds with mass losse s ~ 10~^ M© yr~^ (e.g., 
iHumphrevs fc Davidsonl[l994[ ). Although the line-driven 
winds of pri mordial stars ar c predicted to be weak or 
non-existent (iKrticka fc Kub at 2006), pulsational insta- 
bility of massive stars h as also been fo und to drive mass 
loss (^raffe et al. 2001: iSonoi fc Ume da 2011). Further 
work is necessary to address how massive SMSs could 
form via mass accretion in spite of such disruptive ef- 
fects. 

Stellar evolution under conditions of very rapid mass 
accretion as presented and discussed here is mostly rele- 
vant to the formation of stars in the atomic-cooling halos. 
However, our results could be also important for nor- 
mal Pop HI star formation where H2 molecular cooling 
operates. The typical mass accretion rate for this case 
is around 10~^ Mq yT~^ , but in some exceptional sit- 
uations, e.g., when a progenitor cloud core is extremely 
slow rotating, higher accretion rates M^, ~ 10 ~^ Mq yr~^ 
can be realized (e.g.. IHosokawa et al.l 1201 lb[ ). Since the 
stellar effective temperature is low at ~ 5000 K with 
rapid mass accretion, formation of the HII region would 
be postponed until the mass accretion rate falls below 
10~^ Mq yr~^. This would help the primordial star to 
grow to r nore than 100 M© in th e molecular-cooling halos 
(see also lOmukai fc Pallall2003D . 
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